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A multidisciplinary sensitivity analysis of rotorcraft simulations involving tightly coupled high-fidelity com- 
putational fluid dynamics and comprehensive analysis solvers is presented and evaluated. An unstructured 
sensitivity-enabled Navier-Stokes solver, FUN3D, and a nonlinear flexible multibody dynamics solver, DY- 
MORE, are coupled to predict the aerodynamic loads and structural responses of helicopter rotor blades. A 
discretely-consistent adjoint-based sensitivity analysis available in FUN3D provides sensitivities arising from 
unsteady turbulent flows and unstructured dynamic overset meshes, while a complex-variable approach is used 
to compute DYMORE structural sensitivities with respect to aerodynamic loads. The multidisciplinary sensi- 
tivity analysis is conducted through integrating the sensitivity components from each discipline of the coupled 
system. Numerical results verify accuracy of the FUN3D/DYMORE system by conducting simulations for a 
benchmark rotorcraft test model and comparing solutions with established analyses and experimental data. 
Complex-variable implementation of sensitivity analysis of DY MORE and the coupled FUN3D/DYMORE sys- 
tem is verified by comparing with real-valued analysis and sensitivities. Correctness of adjoint formulations 
for FUN3D/DYMORE interfaces is verified by comparing adjoint-based and complex-variable sensitivities. 
Finally, sensitivities of the lift and drag functions obtained by complex-variable FUN3D/DYMORE simula- 
tions are compared with sensitivities computed by the multidisciplinary sensitivity analysis, which couples 
adjoint-based flow and grid sensitivities of FUN3D and FUN3D/DYMORE interfaces with complex-variable 
sensitivities of DY MORE structural responses. 


I. Introduction 


High-fidelity analysis methods for rotorcraft aeromechanics have become increasingly important to capture com- 
plex physics involving transonic flows, vortical wakes, reversed flows, blade vortex interactions (BVI), and large 
amplitudes of elastic flapping, lagging and torsion motions. A rotorcraft analysis often requires many disciplines such 
as aerodynamics, aeroacoustics, structural dynamics and deformations, flight mechanics, and others to account for 
complex interactions of unsteady fluids with highly flexible rotor blades. Multidisciplinary rotorcraft comprehensive 
analysis (CA) tools! are often used to simulate rotorcraft aeromechanics. These tools are fast and comprehensive, 
but rely on low-fidelity aerodynamic models, such as lifting line and vortex wake models. The state of the art in high- 
fidelity rotorcraft analysis is represented by simulations that couple a CA code with a physics-based, first-principle 
computational fluid dynamics (CFD) code.>~7 

Gradient-based optimization of rotorcraft simulations has been a focus of intensive studies in the past 30 years. -!° 
Most of the studies involving CA tools have been conducted using finite-difference approximations for sensitivities. 
Finite-difference approaches are effective for computing sensitivity of many objective functions with respect to a few 
design parameters. The number of simulations required for such computations is proportional to the number of design 
parameters and may be numbered in hundreds for design and optimization of complex rotorcraft configurations. Finite- 
difference optimization approach may be acceptable for optimization of non-expensive CA models, but it is not feasible 
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for high-fidelity rotorcraft analysis, in which a single large-scale simulation involving a coupled CFD/CA system may 
require hundreds of thousands of CPU hours of high-performance computing. On the other hand, adjoint formulations 
developed in some state of the art CFD solvers enable simultaneous computations of sensitivity of an objective function 
to many design parameters. !'~!4 The computational cost of an adjoint-based analysis approximately equals to the cost 
of a single simulation and does not increase with the number of design parameters. This property makes adjoint 
methods especially suitable for shape optimization where a single or a few objective functions are used, but the 
number of design parameters is large. Discretely consistent adjoint methods offer an additional advantage that the 
computed numerical sensitivities can be rigorously verified through comparisons with the exact sensitivities computed 
by a complex-variable approach. !° In contrast with the finite-difference methods employing real-valued perturbations, 
the complex-variable approach does not suffer from subtractive error cancellation and is capable of producing the 
true sensitivity of the discrete solution with respect to a chosen design parameter. Recently, adjoint-based sensitivity 
analyses of model rotorcraft simulations !!-!®-!7 have been demonstrated to guide design and optimization. 


In the context of multidisciplinary simulations, the sensitivity analysis and the subsequent gradient-based opti- 
mization procedures must include sensitivities across all disciplines, thereby leading to a multidisciplinary sensitivity 
analysis. While discretely-consistent adjoint-based multidisciplinary sensitivity analysis is the ultimate goal, the im- 
plementation of a discretely-consistent formulation in a production-level code is a complex task that may require years 
of development effort. Moreover, development of a multidisciplinary adjoint formulation requires detailed expertise 
and unlimited access to the source codes in each discipline involved. Modern multidisciplinary applications often cou- 
ple large legacy codes developed over many years by different groups and organizations. The codes often exchange 
data through specialized interfaces without exposing the internal treatment of the discipline solutions. In many cases, 
the source codes and/or expertise related to a particular discipline are simply not available. Even if the required ex- 
pertise and source codes are available, the computational time may not be evenly distributed between disciplines. A 
useful multidisciplinary sensitivity analysis tool can couple efficient adjoint-based sensitivities computed for “heavier” 
disciplines such as CFD with sensitivities computed by “black-box” or “grey-box” methods for “lighter” disciplines 
such as CA. The “black-box” methods may involve real or complex-valued finite-difference approaches; the “grey- 
box” methods may involve automatic differentiation methods. !*-!° The sensitivity analysis for each discipline can then 
be integrated within a multidisciplinary sensitivity analysis framework. A theoretical foundation for such integrations 
has been presented in the literature. 7° 


The goal of this paper is to present and evaluate a multidisciplinary sensitivity analysis approach for high-fidelity 
CFD/CA coupled simulations based on a tight coupling methodology. An unstructured, highly-scalable CFD solver, 
FUN3D, 7!” and a nonlinear flexible multibody dynamics CA code, DYMORE,! are coupled to predict the rotor 
airloads and structural responses of helicopter rotor blades. For overset meshes, the domain-connectivity information 
is established via the software libraries described in Ref. 23. A discretely-consistent adjoint-based sensitivity analysis 
available in FUN3D provides sensitivities arising from unsteady turbulent flows and unstructured dynamic overset 
meshes, while a complex-variable approach is used to compute DYMORE structural sensitivities with respect to 
aerodynamic loads. In contrast to the analysis procedure, in which DYMORE operates in the serial execution mode, 
the sensitivity analysis of DY MORE is distributed to all the processors used in the CFD analysis. The computations of 
sensitivity of structural responses with respect to different airload components are embarrassingly parallel and can be 
simultaneously conducted by different processors. The multidisciplinary sensitivity analysis integrates all sensitivity 
components from the coupled systems. 


The material of the paper is presented as follows. Section II highlights the technical approach used in this pa- 
per. A mathematical formulation for adjoint-based sensitivities of FUN3D coupled with complex-variable sensi- 
tivities for DYMORE is presented. Section III provides verification studies for the current implementation of the 
FUN3D/DYMORE (F/D) coupled analysis. The complex-variable analyses for the DY MORE model and the multidis- 
ciplinary F/D model are verified. The adjoint formulations for F/D interfaces and multidisciplinary sensitivity analysis 
are examined and verified by comparing with real-valued finite-difference and complex-variable sensitivities. Finally, 
Section IV summarizes the work and discusses future plans. 


II. Technical Approach 


In this section, the CFD and CA solvers used in this study are first outlined, followed by a brief description 
of the CFD/CA coupling interfaces and the loose and tight coupling procedures. A complex-variable approach for 
computing sensitivities is discussed, and a mathematical formulation of the multidisciplinary CFD/CA sensitivity 
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analysis is given. 


A. Flow and Comprehensive Analysis Solvers 


Solutions of the Reynolds-averaged Navier-Stokes (RANS) equations are computed with the FUN3D flow solver, °?!2? 
developed and supported by NASA Langley Research Center. A standard one-equation Spalart-Allmaras turbulence 
model ”* is used for the current simulations. FUN3D is a finite-volume, node-centered, unstructured-grid RANS solver, 
which is widely used for high-fidelity analysis and adjoint-based design of complex turbulent flows.7>° FUN3D 
solves the governing flow equations on mixed-element grids; the elements can be tetrahedra, pyramids, prisms, and 
hexahedra. At median-dual control-volume faces, the inviscid fluxes are computed using an approximate Riemann 
solver. In the current study, Roe’s flux difference splitting*’ is used. For second-order accuracy, solutions at dual 
faces are obtained by a UMUSCL (Unstructured Monotonic Upstream-centered Scheme for Conservation Laws) 
scheme*®? with the coefficient set to be k = 0.5 for the meanflow equations and « = 0 for the turbulence model. 
The viscous fluxes are discretized with a finite-volume formulation, in which the Green-Gauss theorem is used to 
compute gradients on the dual faces for tetrahedral meshes; this is equivalent to a Galerkin type approximation. For 
nontetrahedral meshes, the Green-Gauss (cell-based) gradients are combined with edge-based gradients to improve 
the h-ellipticity of viscous operators. The diffusion term in the turbulence model is handled in the same fashion as the 
meanflow viscous terms. The vorticity-based source term for the turbulence model is computed using velocity gradi- 
ents evaluated by the Green-Gauss method on dual volumes. To advance the equations in time, a blended second- and 
third-order backward-difference scheme, referred to as BDF2,,°° is employed; the scheme is second-order accurate 
in time, but has a smaller leading error term than the standard second-order backward difference scheme (BDF2). For 
overset meshes, the DiRTlib*! and SUGGAR++7? codes are used to facilitate communications between components 
of the mesh. 

The CA code used in this study is DYMORE, ! a nonlinear flexible multibody dynamics analysis code that provides 
static, dynamic, stability, and trim analyses of rotorcraft configurations. DYMORE contains libraries of primitive 
elements such as rigid bodies, mechanical joints, elastic springs, dampers, beams, shells and plates. For computational 
structural dynamics analysis, DYMORE uses a high-fidelity finite-element method in the time domain without relying 
on a low-fidelity model reduction method. Although DYMORE has not been developed specifically for rotorcraft 
applications, it has been widely used in this area due to its capabilities of handling nonlinear flexible systems with 
arbitrary topologies. The internal aerodynamic model in DYMORE is a low fidelity approximation based on lifting 
line theory and vortex wake models. To enable sensitivity computations conducted via a complex-variable approach, 
a complex-variable version of DY MORE has been developed and verified by comparing with real-valued analyses. 


B. CFD/CA Coupling Approaches 


An F/D interface has been previously developed to enable the exchange of aerodynamic loads and structural deflections 
between high-fidelity aero and structural dynamics models. ** In the current work, the interface has been extended to 
handle mixed-mode communications between real-mode FUN3D and complex-mode DYMORE as well as to enable 
multidisciplinary sensitivity analysis. The new implementation of the F/D system can conduct loose or tight coupling 
simulations. 

In a loose coupling approach, *? an initial solution from DYMORE, using its internal linear aerodynamics model, 
provides an initial estimate of the elastic blade deformations for FUN3D. CFD computations start from a freestream 
flow; delta airloads between the CFD solution and lifting-line computation in DYMORE are collected along with the 
time advancement process for one revolution in the first coupling iteration and 2/N, revolutions in the subsequent 
coupling iterations. Here, N, denotes the number of rotor blades. The use of one full revolution in the first coupling 
iteration aims to reduce the transient effects caused by freestream initial conditions. Once the computation of the flow 
field is complete, the delta airloads are passed to DY MORE to compute trimmed elastic motions, which are transferred 
back to FUN3D to generate new blade surface geometries and motions. This process is repeated until a converged 
periodic solution is achieved. 

In a tight coupling approach, the CFD and CA models are still executed separately, but the data exchange occurs 
once per time step. Fig. 1 shows a diagram illustrating the F/D tight coupling procedures, in which FUN3D functions 
that have DYMORE structural solutions as inputs are highlighted in red. At each coupling iteration, DY MORE 
performs one time-step advancement with a set of sectional airloads as inputs. The structural solutions describing 
linear and angular displacements at the quarter-chord line of rotor blades are sent to FUN3D. High-order surface 
spline functions are created, followed by an extraction of average rigid motions from the CA displacements and a 
reposition of the CFD surface grid. The CFD volume grid is rigidly moved and elastically deformed to accommodate 
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Figure 1. Diagram illustrating the tight coupling procedure of F/D models. 


the new surface mesh. The CFD flow problem is solved on the updated mesh. When a desired convergence level 
is reached, a boundary slicing procedure is performed to obtain updated sectional airloads (forces and moments) at 
specified airstation locations. The new CFD airloads are then transferred to DY MORE through the interface to start 
the next coupling iteration. In the present implementation of the tight coupling procedures, the data exchange between 
the CFD and CA models is performed in memory by linking FUN3D with a prebuilt DY MORE library, which results 
in direct, fast data access. For verification purposes, the tight coupling computation is conducted with trim control 
angles that have been obtained from a loose coupling procedure. The trim control parameters are not updated during 
the tight coupling analysis procedure, and therefore trim computations are avoided. This simplification significantly 
accelerates the iterative convergence of DYMORE solutions. In future design and optimization applications, the trim 
control parameters will be updated through constraint optimization procedures. !°?° In Section III.A, a verification 
study is conducted to demonstrate the accuracy of the current F/D loose and tight coupling procedures for the HART- 
II baseline rotorcraft configuration. ’ 


C. Complex-Variable Approach 


In this work, the complex-variable method !> is employed by the multidisciplinary CFD/CA sensitivity analysis to 
compute sensitivities of structural responses with respect to aerodynamic loads and to verify multidisciplinary sensi- 
tivities. A brief description of the method is given as follows. 

For a real-valued function, f(x), of a real variable x, the complex-variable method allows simultaneous evaluation 
of the function and its derivative, df /dx. For an input, x+ ih, with an imaginary perturbation, ih, the complex-valued 
function can be expanded in a Taylor series as: 


df Wd f eat 


ih) = + ih ae 1 
f(xtih) =f) +ihe — 7 aa te Ga (1) 
From this equation, the real part of the complex-valued expansion gives the function itself, 

F(x) © Re[f(x+ ih)], (2) 
and the imaginary part yields the derivative 

df Iml|f(x+ih)] (3) 

dx h d 
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with a truncation error O(h7). 

Unlike a finite-difference method with a real-valued perturbation, the complex-variable approach does not suffer 
from subtractive cancellation errors, thereby providing increased accuracy and robustness. Within machine accuracy, 
derivatives computed by the complex-variable method with perturbation h = 10~>° can be considered as the true 
derivatives with zero truncation error. 

The standard FUN3D software suite allows for complex-mode operation to verify hand-differentiated sensitivities, 
but a similar capability was not available in DYMORE. For the present study, a complex mode has been developed in 
DYMORE by replacing all double and float variables declared in the original (i.e., real-valued) DYMORE code with 
complex variables of double precision and redefining operators needed for complex-valued operations. The complex 
mode of DYMORE enables the evaluation of comprehensive analysis sensitivities via the complex-variable approach. 
The F/D interface has been extended to operate in the complex mode, supporting the exchange of complex-valued 
aerodynamic forces, moments and surface displacements. 


D. Sensitivity Analysis 


A mathematical formulation for the sensitivities of a tight coupling multidisciplinary CFD/CA analysis is given in this 
section. The following notations are used: u” is a CA solution vector at a given time level n. The vector is used to 
update the surface position and the average rigid motion of component CFD grids at time level n. Initial CA solution, 
w°, is given as a solution of a loose coupling formulation and is independent of design parameters. A vector of airloads, 
f”, is computed by the CFD code at time level n at specified airstations. The fixed baseline volume and surface grids, 
X and X,, are generated based on the baseline blade configuration, independent of design parameters, and used to 


specify the baseline elasticity matrix, K, 


K X= X,. (4) 


The reference component volume and surface grids, X and X,, account for variations of shape design parameters (used 
to perturb the surface geometry), D, and define the reference frame. The reference blade grid X is computed by solving 
a mesh elasticity problem 


KX = X,(D). (5) 


The overset volume and surface grids, X” and X”%, reflect the actual positions of all blades at time level n. The initial 
grid, X°, corresponding to time level n = 0 is computed by rigidly moving the reference blade grid to the corresponding 
initial position for each blade and combining these component grids with a stationary background grid. 


X°=7° X, (6) 
po _ PS 70) Ri0) G/q,0 
T, =T (8(u’), B(u"), O(u"))To,. (7) 


Here, 5,8, and 0 are lead-lag, flap, and pitch control angles, respectively; 5, B, and 6 denote quantities averaged over the 
entire span of a rotor blade; T’(5(u°), B(u°) ,6(u°)) is an extracted (composite) rigid motion corresponding to averaged 
CA displacements. As described in Ref. 32, this averaged rigid motion helps keep the elastically-deflected blade 
centered within the blade volume mesh, and helps prevent negative volumes during the mesh deformation process, 
when large flapping motions are present. The transform matrix for the composite rigid motion is obtained by successive 
multiplications of the transform matrices corresponding to each rigid motion, 


T(5,B,6) = T(5)- T(B) -T(8). (8) 
Each rigid motion is represented as a 4 x 4 matrix 


Ry Ro Rp tk 
Ro, Ro R23 ty 
R31 R3o R33 it 

0 0 0 1 


(9) 


Here, the 3 x 3 block matrix R defines a general rotation and the vector t appearing in the fourth column of the 
transform matrix defines a translation. The motion Tg, in Eq. (7) is a rigid motion transforming a blade from its initial 
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position to the reference frame, in which the CA displacements are defined. The grid equations at times n > 1 include 
elastic grid deformations satisfying 


A 


K"X" = T" X,(D), (10) 
Ty = Ty T"(8(u"), B(u"), @(u") To, (11) 


Here, T, denotes the elastic motions applied to individual surface nodes, and Ty is a rigid motion corresponding to the 
azimuthal rotation of the blades. The elasticity matrix, K”, is derived from a rigidly moved reference grid, satisfying: 


i Xi (12) 
wT (8(u"),B(u"), 6(u")) To, (13) 


Note that the elasticity matrix at time level n relates to the baseline elasticity matrix as 
nap Ana 
KT =T'K. (14) 


The vector Q” is a vector of CFD solutions at time level n. The solution at time n = 0 is specified from a precomputed 
loose coupling solution, Q, that may depend on CFD design parameters such as angle of attack, Mach number, etc. 
Given variables that are independent of design parameters, K, u°, and To), and variables that may depend on design 
parameters, Q(D), X,(D), and Ty(D), the following equations for a tight coupling formulation are written in the order 
in which the equations are solved: 


Reference Grid: G(X,X,,D) =K X— X,(D) =0, (15) 
Initial Grid: G°(X°,X,T?(u°),D) =X°_T° X=0, (16) 
Initial Flow: R°(Q°,Q(D)) =Q°— Q(D) =0, (17) 
Initial Load: F°(f°, Q°,X°) =0, (18) 

CA: C'(u",u!/£""!D) =0, (l<n<N) (19) 

CFD Grid: G"(X",X;,T)(u"),T"(u"”),D) — =T) K(T")-'X" —T’X,(D) =0,  (I<n<N)_ (20) 
CED Flow: R”(Q”,Q”!,x",x”"! D) =0, (l<n<N) (21) 
Loads: F"(f"”,Q”,X”) =f" — §"(Q",X") =0. (l<n<N) (22) 


Here, S denotes a boundary slicing function performed in FUN3D for obtaining sectional airloads. For simplicity of 
presentation, the time-dependent equations are assumed to depend on solutions at the current and one preceding time 
level only. In the actual solver, the time derivatives are discretized through second-order backward difference schemes 
that use solutions at the current and up to three preceding time levels. The residuals of the CFD equations are in the 
form of an Arbitrary Eulerian-Lagrangian formulation that is suitable for dynamic deformable grids and accounts for 
grid speeds and geometric conservation law. The details of CFD equations on overset grids are provided in Ref. 25. 

In the sensitivity analysis framework presented here, the CFD flow and grid equations and the loads equations use 
adjoint formulations, whereas sensitivities of the surface-grid elastic motion T, and the extracted rigid motion T, to the 
aerodynamic loads f” are computed by the complex-variable method. For a general objective function g(Q,u, X,f,D), 
the adjoint formulation is derived via the Lagrangian functional formed as 


N N N 
L=g(Q.u,X,f,D) +) (Ag)7R"+ Y (AY) G" + YAR F" + [Ag] (23) 
n=0 n=0 neal 

where Q,u, X, and f represent CFD solutions, CA solutions, grid, and load solutions at all time steps, respectively; Ar, 
AG, and Ar are the vectors of the time-dependent adjoint solutions for the flow, grid and loads equations, respectively; 
Ag denotes the adjoint solution for the reference grid; and the superscript T represents the transposition operator. 
Differentiating the Lagrangian with respect to design variables D and equating the coefficients of state variable sensi- 
tivities to zero** yields the following adjoint equations: 
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Grid : 


Load : 


Fa + [AF] Ser =0, 
2s] + on” [2p] + oer [2a] + an” [2 -», 
ae] +oer Be] 3, [el ( e mr, aps S|) -0, 
oe] + tan? [SE] + taney” [SE] + tani” |] ~0, 
2 | + tact” |S] + tan? [S| + tage)” [| san [SE] =o, 
ak 
lms] Sia (](F-Celee) 
Ea + [Ag] ES + [AR] so + [Ap] san =0, 
i] Son oon] 
inal” |e] + B98” [2] 0 


(24) 


(25) 


(26) 


(27) 


(28) 


(29) 


(30) 


(31) 


(32) 


(33) 


The adjoint equations are solved in the inverse order in time. In the current implementation, [dF” /df”] is the iden- 
tity matrix, thus the load adjoint equation becomes a simple algebraic equation. The summation terms appearing in 
Egs. (27) and (30) correspond to the sensitivities of rigid motions and surface elastic motions to sectional airload vari- 
ations, which are computed via imaginary perturbations of individual components of the airloads f” at each time level 


n. For each such perturbation, the complex-variable CA solution and the transform matrices, T and 1 are computed 
for time levels k,n <k <N. Incontrast to the analysis procedure, in which the CA code operates on a single processor, 
the computations of sensitivities to airload perturbations are distributed to all the processors used in the CFD compu- 
tation. The total computational time for evaluating the complex-variable sensitivities is proportional to the number of 
airstations that control the data exchange between the structural and aerodynamic models, the number of time steps 
N, and inversely proportional to the number of processors used in parallel computing. The contributions of struc- 
tural sensitivities to the CFD flow and grid adjoint solutions are accounted for through the vectors [A”.]’ [OF”/0Q”| 
and [A”.]’ [9F”/dX"] in the flow and grid adjoint equations, respectively. Here, [F” /dQ"] and [dF"/dX”| represent 
the sectional airloads sensitivities with respect to the conserved flow variables and grid solutions, respectively. A 
verification test for the sectional airloads sensitivities is provided in Section III.C. At time level n = 0, the matrices, 
[OF° /of?], [OR°/dQ°], and [dG°/dX°], are identity matrices, and the load, flow and grid adjoint equations become 


simple algebraic equations. 
The general form of the sensitivity derivatives of the objective function is written as: 


dL _ [ag 0, [8G°] [00] Xr fOR" ot { [aG°] [dX,]  [aG°] aT? 
i * Se] + (asl BE ap|* 2,44) lap] * el | lax] [ap] * [ag] ao] ) SO 
s dG") [dX dG"] | oT" dG" ] oT” r (| 0G | [dX 
nj Sly 7 Ss a S 
+ Sioa" ((5e] [Se] + [5e] [Sn] + [Bt] [ao] +a" (| [0] 
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Specifically, the sensitivity with respect to flow design parameters (e.g., angle of attack) is computed as 
dL dg 06°] [oQ]) JX 7 [OR" 
a Apg| | ==] | =~ Ar) | =~]; 
dD lap] ian Eq Es + DA Lap | er 


the sensitivity with respect to shape design parameters is computed as 


A 


0G 
OX; 


dL dG" | | dX; T 
ag x —1 + |Aa 

dD =| >| a al [Se | [a Ae 

the sensitivity with respect to kinematics design parameters (e.g., rotation rate) is computed as 


Zea (Ge] [35]+ Ge] Ge])- 


n=1 


dX.) 
4). 06 


aT” 
aD 


ag? 
aT? 


ar" 
oD 


dL = og 017 
7D > B + [AG] 


Ill. Numerical Results 


This section provides verification of the current implementation of loose and tight coupling procedures for CFD/CA 
systems by comparing the solutions for a test model of the baseline HART-II rotorcraft configuration’*> with estab- 
lished benchmark computations and experimental data. The complex-variable analysis is verified by comparisons 
with relevant real-valued, finite-difference computations. The adjoint formulation for sectional airloads computation 
is verified with the complex-variable approach. Sensitivities of the complete multidisciplinary CFD/CA computations 
for lift and drag functions to the angle-of-attack design parameter are presented. 


A. HART-II Rotorcraft Analysis 
I. CFD/CA Loose Coupling Simulations 


The flow conditions in the HART-II baseline case include a tip Mach number of 0.6387, an advance ratio of 0.151, 
and a shaft tilt angle of 4.5° (adjusted after the wind tunnel wall correction). This case represents a descending flight, 
where blade-vortex interactions are significant. Other important variables for the blade geometry and test conditions 
are detailed elsewhere.’ Solutions for the HART-II baseline case were obtained with an unstructured mixed-element 
mesh containing approximately 14 million nodes and are compared to established solutions of FUN3D/CAMRADII 
(NL-2)° and a previous implementation of FUN3D/DYMORE (GIT-1)’ as well as experimental data. In the present 
study, the BDF2,,, temporal scheme is used to advance the flow computation in time with a time step corresponding 
to the increment of 1° azimuth angle. In ten coupling cycles that involve 5.5 revolutions, the variations in the control 
angles and thrust deltas are within 0.01° and 0.1%, respectively. Figure 2 illustrates vortex structures near the HART- 
II rotor in the baseline case, where the BVI phenomena and the formation of blade tip vortices are clearly observed. 
Converged solutions of aerodynamic loading and blade tip motions are examined in the following. 

Figures 3(a) and (b) demonstrate comparisons of the normal force and pitching moment at the 87% blade radial 
station evaluated from the current approach (F/D), NL-2, and GIT-1 solutions, and the measurement data. The mag- 
nitudes and shapes of the normal force predicted by the current F/D solution are in good agreement with those of the 
NL-2 solutions. Furthermore, the events of BVI at the advancing side and particularly the retreating side are cap- 
tured well by the F/D analysis. The under prediction between 5° and 40° and the over prediction between 230° and 
270° appearing in the GIT-1 solution are not present in the current F/D simulation. The differences observed here are 
possibly due to the fact that the present implementation of F/D coupling procedures follows closely with those of the 
FUN3D/CAMRADII development to date, while the GIT-1 solution was obtained with an earlier version of FUN3D, 
where differences exist in the global slicing procedures for extracting the blade sectional airloads and the conversions 
of elastic motions to the blade surface geometries. A similar comparison is made for the aerodynamic pitching mo- 
ment, where the solutions of F/D, NL-2 and GIT-1 methods are generally in a good agreement, particularly in the 
retreating side. 

Figure 4 compares the azimuthal variation of elastic motions including flap, lead-lag and torsion deformations at 
the blade tip. The flap deflections are obtained by removing the pre-cone angles (2.5°) from the vertical displace- 
ments, while the elastic torsion is obtained by excluding the pre-twist angles and the pitch control inputs from the 
total geometric pitch angle. As demonstrated in Fig. 4(a), the overall characteristics of the predicted blade tip flap 
deflections in the F/D solution agree well with those of the NL-2 solution, both in magnitude and phase. Compared 
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Figure 3. Comparisons of normal force and pitching moment predictions at 87% radial station with means removed for the baseline 
HART-II rotor. 


to the GIT-1 solution, some differences can be observed, especially in the first and third quadrants of the rotor disk. 
These discrepancies are likely caused by the differences in the airload predictions as observed in Figs. 3(a)-(b). All 
the computational results exhibit a secondary elastic flap maximum at approximately 40°-50° azimuth, while the 
magnitude of the current solution is higher than that of the GIT-1 result but quite similar to that of the NL-2 solution. 
Blade tip lead-lag deflections are evaluated in Fig. 4(b), where consistent agreement between motions of the F/D and 
NL-2 solutions is again achieved. Compared to the measurement data, a constant offset roughly amounting to 1/5 of 
the chords in the elastic lead-lag deflections is also observed, but this translational offset behavior is commonly seen 
in other work.*°77 Figure 4(c) compares the elastic torsion at the blade tip, where the 2/rev behavior of the elastic 
twist is predicted well in the F/D solution. In addition, the overall torsion levels are more accurately captured by the 
F/D solution as compared to the GIT-1 solution, which shows approximately 1° under predictions across all azimuth 
locations. 


2. CFD/CA Tight Coupling Simulations 


The tight coupling F/D algorithm is verified in this section by comparing tight and loose coupling simulations con- 
ducted with the same trim control angles. Trim computations are generally performed with a loose coupling CFD/CA 
approach. To achieve converged trim solutions, the CA model typically simulates many revolutions per each coupling 
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Figure 4. Predictions of blade tip flap, lead-lag and torsion deflections for the baseline HART-II rotor. 


iteration, which becomes expensive for a tight coupling procedure. In the present work, the trim control angles are 
taken from a fully converged loose coupling solution. The tight coupling simulations then start from the loose coupling 
solution, and the collective and cyclic pitch control angles prescribed by the loose coupling trimmed solutions are kept 
throughout the tight coupling computation. Such a tight coupling procedure should result in minimal change in the 
solutions. 

Figure 5 shows comparisons between the loose and tight coupling F/D solutions for predictions of the normal force 
and pitch moment at 87% of the radial station, and the tip flap deflections in the HART-II baseline rotorcraft test case. 
One can observe that the tight and loose coupling procedures deliver almost identical solutions, which implies that the 
tight coupling F/D solver functions correctly, and the rotor maintains the trimmed condition. 


B. Complex-Valued CFD/CA Analysis 


Complex-variable sensitivities for the coupled CFD/CA analysis can be obtained by operating both models in fully 
complex modes and perturbing the imaginary part of a design parameter. One such computation can be used to obtain 
the sensitivities of multiple objective functions with respect to a single design parameter. To ensure that the complex- 
valued DYMORE model functions properly, real parts of the complex-valued DY MORE solutions describing the linear 
and angular displacements at the quarter-chord line of a rotor blade in the HART-II baseline rotorcraft test case are 
examined. The examination shows that the real parts of the DY MORE solutions are identical (matching all significant 
digits) to those obtained by the real-valued DYMORE model. 

Table | lists aerodynamic lift at time level n = 3, computed by the real-valued F/D tight coupling procedures and by 
the complex-valued F/D tight coupling procedures with several perturbed design parameters including angle of attack 
and shape design parameters on different blade surfaces of the HART-II baseline rotorcraft configuration. The real 
parts of the lift function in the complex analyses agree with that computed by the real-valued F/D analysis to machine 
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Figure 5. Comparisons of normal force and pitching moment predictions at 87% radial station with means removed for the baseline 
HART-II rotor with tight and loose coupling approaches on a 7M node unstructured mesh. 


Table 1. Real-valued and complex F/D analyses for lift (at time level 3) in the baseline HART-II rotorcraft test case. 


Analysis __Perturbed parameters _— Real part of solution Sensitivity 
Real 0.0007989 18217192 
Angle of attack 0.0007989 18217192 0.000024054292434 


Shape variable, blade 1 0.000798918217192 —0.000044901491848 
Complex Shape variable, blade 2 0.000798918217192 0.000009593294576 
Shape variable, blade 3. 0.000798918217192 | —0.000000293789365 
Shape variable, blade 4 0.000798918217192 0.0000058 15724362 


zero. The imaginary parts of the objective functions in the complex analyses correspond to the lift sensitivities with 
respect to the specified design parameters and can be used for verification of the multidisciplinary CFD/CA sensitivity 
analysis. 
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Table 2. Comparisons of sectional airloads sensitivities with respect to the grid coordinates (y-component) and conserved flow variables 


(x-component of momentum) at 70% radial station via direct differentiation and complex-variable approaches. 


Method oF,./dQ dF, /dQ OF, /dQ 
Adjoint 0.000041123560446 0.000126949885201  -0.006184933414864 
Complex variable 0.000041123560446 0.000126949885201 -0.006184933414864 
Methods OF inx/dQ OF ny /OQ OF nz/OQ 
Adjoint -0.000035570041307 -0.000005662910700 = -0.000000350458662 
Complex variable -0.000035570041307 -0.000005662910700 = -0.000000350458662 
Methods OF, /0X OF, /0X OF, /dX 
Adjoint -1.66200987502660 1.97008265283363 0.520855236517486 
Complex variable = -1.66200987502690 1.97008265283398 = 0.520855236517574 
Methods OF nx /OX OF ny /OX OF nz /OX 
Adjoint -0.015116700097615  -0.015298667011763 0.007316138436884 
Complex variable -0.015116700097618 -0.015298667011765 0.007316138436883 


Table 3. Sensitivities of lift and drag to the angle-of-attack design parameter evaluated by the real-valued central difference (CD) meth- 
ods with perturbation sizes of 10~° and 10~°, complex-variable approach with perturbation size of 10~>°, and the coupled F/D adjoint 


sensitivity analysis approach. 


Approach Lift sensitivity Drag sensitivity 
Real-valued CD, h = 10~° 0.000024054294350 0.000015909433903 
Real-valued CD, h = 10~° 0.000024047540541 0.000015921931245 

Complex variable, h = 10~°° 0.000024054292434 0.000015909425336 
Multidisciplinary adjoint 0.000024054286300 0.000015909424183 


C. Coupled CFD/CA Sensitivity Computation 


The complex-variable method is employed to compute sensitivities of structural responses with respect to aerodynamic 
loads and relevant design parameters. All other derivatives required for the evaluation of the adjoint variables and 
the final sensitivity are implemented by hand differentiation of the corresponding routines in the CFD/CA system. 
In this section, a verification test for sectional airloads sensitivities with respect to shape variables and conserved 
flow variables is performed via the complex-variable approach for the HART-II rotorcraft test example discussed in 
Section III.A. Table 2 shows comparisons of representative aerodynamic airloads derivatives at 70% radial station 
computed by the hand differentiation and complex-variable approaches for a given state of viscous flow and grid 
solutions. As observed, the hand-differentiated sensitivities agree with complex-variable sensitivities well; at least 
12-digit agreement is obtained, which is considered as a sufficiently good agreement.”° 

The sensitivities of aerodynamic lift and drag functions evaluated at time level n = 3 with respect to the angle 
of attack is computed by the multidisciplinary F/D sensitivity analysis that integrates adjoint-based sensitivities of 
FUN3D and F/D interfaces with complex-variable sensitivities of DY MORE (cf. Eq. (35)). Although only the flow 
adjoint solutions are required in the evaluation of the sensitivity to a flow design parameter, the grid adjoint solutions 
must also be computed due to the coupling. For verification, all the time-dependent problems in both forward and 
adjoint analyses are converged to machine zero residuals. For demonstration purpose, the first-order backward differ- 
ence temporal scheme is employed in this example; other higher-order temporal schemes available in FUN3D will be 
tested in the future work. 

Table 3 lists the sensitivities of lift and drag at time level n = 3 to the angle-of-attack design parameter computed by 
the real-valued finite difference (second-order central difference) methods with perturbation sizes of 10~° and 10°, 
respectively, the complex-variable approach with a perturbation size of 10~>°, and the multidisciplinary F/D sensitiv- 
ity analysis. Compared to the complex-variable sensitivity, the real-valued finite-difference results corresponding to 
the perturbation size of 10~° show good accuracy in the estimations of lift and drag sensitivities, while the accuracy 
degrades significantly as the smaller perturbation size is selected. This loss of accuracy of the finite-difference deriva- 
tives with small perturbation size is expected due to subtractive cancellation errors, and highlights the shortcomings of 
real-valued finite differences. The multidisciplinary F/D sensitivities show good agreement with the complex-variable 
results; approximately 11 digit match is achieved in both lift and drag sensitivities. 
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IV. Conclusions 


This paper has presented a multidisciplinary sensitivity analysis framework for high-fidelity rotorcraft computa- 
tions conducted by a tightly coupled multidisciplinary system including a computational fluid dynamics code, FUN3D, 
and a rotorcraft comprehensive analysis code, DYMORE. The implementation of the current FUN3D/DYMORE 
(F/D) analysis procedures has been verified by simulating the baseline HART-II test model. The computed solutions 
show a good agreement with established predictions of unsteady aerodynamic airloads and various elastic deflections. 
Discretely-consistent adjoint formulations have been developed and verified for F/D interfaces. Structure-related sen- 
sitivities of deforming and rigid motions of component grids are evaluated through a complex-variable approach. The 
complex-variable method verifies sensitivity derivatives computed with the adjoint-based methods. Complex-variable 
implementation of sensitivity analysis of DYMORE and the coupled F/D system is verified by comparing with real- 
valued analysis and sensitivities. Results show that the blade sectional airload derivatives computed by the adjoint- 
based and complex-variable methods are in good agreement. Sensitivities of the lift and drag functions obtained by 
complex-variable F/D simulations are compared with sensitivities computed by the multidisciplinary sensitivity anal- 
ysis, which couples adjoint-based flow and grid sensitivities of FUN3D and F/D interfaces with complex-variable 
sensitivities of DYMORE structural responses. Future work will be concentrated on the verification of the coupled 
F/D multidisciplinary sensitivity analysis approach for all types of design parameters and applying it to gradient-based 
multidisciplinary optimization of rotorcraft configurations. 
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